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I. INTRODUCTION 



A. Motivation 

It is remarkable how little is known about classical solutions of relativistic many-body 
electrodynamics with radiative damping. Textbooks often describe simple non-relativistic 
solutions to the one-body problem with an external force and radiative damping. Studies in 
accelerators and antennas also consider radiation emitted from charged particles. However, 
when many bodies are included and the back-reaction of the radiation on the particles 
involved is considered, the resulting equations are amenable neither to analytic solutions 
nor to numerics. We explain the reason for this in greater detail below but we emphasize 
here that the difficulties involved are not merely technical but relate to important questions 
about the nature of classical electrodynamics. 

In fact, until recently, there has been no systematic development of techniques to integrate 
the equations of relativistic many-body electrodynamics with radiative damping. One reason 
for this is that, hitherto, there has been no compelling physical reason to consider the 
intricacies of this problem. For example, in astrophysical processes, the radiation emitted 
by charged bodies is of interest; however the back-reaction of this radiation on the bodies is 
small compared to other forces. Hence, the radiation reaction is often ignored or included 
only heuristically. 

In recent years, this situation has changed. Simulations of gamma ray bursts, for example, 
have reached a level of sophistication that authors have started to investigate the influence 

n 

of radiation reaction 111. A more immediate motivation for this study comes from papers 

□ n n 

of Umm- In these 

papers, it was conjectured that some classical orbits of the many-body 
system have special properties when radiation reaction is considered. However, further 
progress on this front was prevented by the absence of techniques to integrate the many- 
body system 1 with radiation damping. 

In this paper, we develop such a technique. Our technique may be applied to check the 
conjectures of Moreover, this technique is of independent interest because it is 

applicable to the investigation of radiative damping effects in astrophysical processes. It 
would also be interesting to include the effect of radiative damping in molecular-dynamic 



1 We explain the reason for our emphasis on many-body in Section [I CI below 
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simulations carried out in codes such as charmm fl, fl, AMBER(7| or GROMOS j 

In this paper, we elaborate our technique and explain its advantages. Furthermore, we 
bench-mark our technique by applying it to the one-body problem. This is feasible since, 
as we explain further in section ITCl some solutions to the one-body Lorentz-Dirac equation 
have already been obtained. Our technique agrees, excellently, with these solutions. Having 
established the validity of our methods, we hope to apply them to the concrete physical 
problems described above in a subsequent paper. 



B. Setting 

Consider a point particle with mass m, charge q and world-line a^{r) parameterized by 
its proper time r. In the presence of an external electromagnetic field strength Fj^ the 
equation of motion of this particle is given by: 

J r wz, i 
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The second term on the right may be understood as the regularized action of the particle 
on itself. We emphasize that the structure of this term is quite delicate and it cannot be 
modified without modifying the Maxwell theory itself. 

This is a third order differential equation and requires 3 initial conditions. However, for 
a generic choice of initial conditions a(0), d(0), a(0), we obtain a runaway solution where 
the velocity of the particle increases continuously in the future! This solution is clearly 
unphysical and the existence of such runaways is an embarrassment for the classical theory 
of electromagnetism. Several solutions to this problem have been proposed in the literature. 
We discuss some of these, in Section II Dl and Appendix [S] However, one common solution 
to this problem is as follows. Let us say, we assign the initial acceleration to satisfy 



«a(0) 
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with k = ■^- 2 — 5-, then the asymptotic solution is free of runaways 91. 

A commonly articulated philosophical problem with this is that it requires one to know 
all forces that act in the future on the particle. Moreover, the particle starts accelerating 
before any force is actually applied. This is the problem of pre-acceleration. 

However, we are more concerned with a different problem that makes Formula (|2J) com- 
pletely unsuited to any actual computations. The point is this: if we choose a(0), to be even 



infinitesimally different from the value given by (J2]), the solution ceases to be bounded and 
turns into a violent runaway. Stated another way, equation (j2]) specifies a co-dimension 1 
surface in the space of all possible initial conditions. The behavior of solutions to equation 
(CE]) is discontinuous in this space. On this surface, the solutions are bounded; away from 
this surface, the solutions comprise runaways. 

From a computational point of view, this makes any numerical computations with (J2]) 
impossible since we need to specify initial conditions exactly at the critical value prescribed 
by (J2]). From a philosophical point of view, this need to 'fine tune' the initial conditions to 
obtain a meaningful solution is rather unappealing. 

Hence, in this paper we do not adopt this point of view but look for another method to 
tame the runaway solutions of classical electrodynamics. 



C. One-Body vs Many-Body 

We now pause to explain a critical difference between the many-body system and the 
one-body system. For the one-body system there exists a simple method that often allows 
us to sidestep the above complications. 

This involves solving equation ([1]) backward in time. Here, the runaways solutions get 
rapidly damped, leaving behind the 'physical' solution. To our knowledge this technique 
was first used by Plass to obtain the solution for a charged particle orbiting an infinitely 
massive object with the opposite charge [lo| . 

However, for the many-body problem this method fails. This is because the field acting 
on a particle at a given time depends on the positions of the other particles in the past. 
Stated another way, the equations of motion for a relativistic many-body system are func- 
tional differential equations. Such equations can be solved forward in time, given past data. 
However, in general, they do not have a unique solution backward in time. 2 Hence the trick 
of solving the equations of motion backward in time to get rid of the runaways cannot be 
applied to the many-body problem. 

In this paper, we will develop a technique that is applicable to the many-body system. 
As we have already mentioned, we will check this technique against extant solutions of the 



2 This raises some difficulties for the method discussed in 
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one-body problem obtained using the trick above. Thus this trick serves to bench-mark our 
method; while the trick itself does not generalize to the many-body case, our method does. 



D. Other approaches to this problem 

The formidable conceptual and calculational problems outlined above are sometimes 
wished away by appealing to quantum mechanics. This argument goes as follows. Equa- 
tion (TjQ) is only an approximate equation valid at distances large compared to the Compton 
wavelength of the particle. At shorter distances we should use the full quantum theory. 
Hence, according to some authors, we need not take the difficulties above seriously, since 
the equation that they arise from is only approximate in the first place. 

We do not find this point of view satisfactory. For one, we believe that the classical 
theory should be well defined in its own right. Second, in the concrete physical situations 
where we wish to use (JTJ) , it is hardly feasible to use the entire quantum theory. 

However, the observation that equation (TjQ) is valid only at distances large compared to 
the Compton wavelength suggests another possible resolution to the problem. This is to 
modify itself at distances that are short compared to the Compton wavelength of the 
particle. 

This approach has been adopted by many authors. One popular idea has focused on 



giving the particle itsel 
we refer the reader to 



a structure. For an excellent review of models of extended particles 



121 ] . However, it is hard to construct a Lorentz invariant theory of 



13|). one needs to worry about 



extended particles. When this can be done (see for example 
the forces that keep the extended particle stable. In fact, additional terms coming from 
the energy due to the internal stresses that hold the particle together, must be added to 
the mass of t he p article to obtain covariant equations of motion. This is the famous '4/3' 
problem (see [131] for a review). 

A second approach is to place a momentum cutoff in the theory. Since a momentum 
cutoff can be only be placed after Wick rotating the momenta to Euclidean momenta, it 
respects neither Lorentz invariance nor gauge invariance. 

We discuss a third approach to this problem in Appendix [A] 
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E. Our Technique 



To summarize the discussion above, we are looking for a modification of equation ([I]). 
We would like this equation to be a natural consequence of a modified classical theory. 
Furthermore, we would like such a theory to satisfy four properties: 

1. It should be Lorentz covariant and gauge invariant. 

2. Given generic past data, the theory should possess regular solutions in the future. 

3. It should reduce to the usual Maxwell theory at large distances 

4. It should be simple enough to make solutions of the many-body problem with radiation 
damping feasible. 

We have seen that the method of assigning the initial acceleration via equation ([2]) fails 
on count [2] above while other methods of assigning a structure to the particle fail on counts 
[TJ or |H Yet other methods discussed in Appendix |A] are based on ad hoc approximations 
that are unacceptable because they modify the second term in ([T]) without modifying the 
Maxwell theory. In this paper, we will propose a modification of Maxwell electrodynamics 
that satisfies all four conditions above. 

We modify the propagators of classical electrodynamics at small scales in position space, 
while retaining Lorentz covariance and gauge invariance. The effect is as if we had a 
covariant model of extended particles. Since, however, it is electrodynamics which is being 
modified, this allows us to sidestep issues of internal stability of the particle and the equations 
of motion are covariant without the addition of ad hoc terms. Further, the simplicity of this 
model makes numerical calculations feasible. At long distances our theory reduces to the 
usual Maxwell theory. 

In the analogy with extended particles, the length-scale at which electrodynamics is mod- 
ified corresponds intuitively to the 'size' of the extended particle. This finite size introduces 
a delay in the interaction of the particle with itself and so we obtain a functional differential 
equation even for the one-body problem with radiation damping. Although the resulting 
FDEs are numerically stiff, we are able to solve them using numerical techniques recently 
developed by Guglielmi and Hairer 

w. 
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We find the following surprising result. Locally, the force of radiation damping is well 
approximated by the customary third order expression. However, globally our FDE has 
dramatically different properties from the solutions of equation ([IT) . In the appropriate 
regime of parameters, our theory does not exhibit any runaway behavior. The first part 
of our result ensures that standard calculations of radiation emitted by antennae are not 
affected noticeably. 

As already mentioned, the simplicity of this theory makes numerical solutions of the 
many-body problem feasible. In this paper, we solve the important test case of a radiating 
particle in a central potential. The particle gradually spirals in, as is expected on physical 
grounds. Furthermore, we are able to read off the rate of the decay of this orbit. As 
already stated for this one-body problem, we can obtain this solution, even in Maxwell 
electrodynamics, by means of the trick explained in Section II CI The two solutions are in 
good agreement. 

In a subsequent paper we hope to apply the techniques developed here to the full 
relativistic two-body problem. 



II. THE MODEL 

Given a 4-current density, j^(x), the potentials, are given by the formula 

V x ) = ^ / ^(y)G(x, y)rf 4 y + d, X , (3) 

where the (retarded) Green function is given by [ljj 3 

G(x, y) = 8(\x — y| 2 )6'(x — y°). (Maxwell electrodynamics.) (4) 

X is an arbitrary scalar function which we will henceforth ignore by choosing Lorenz gauge 
d^A^ = 0. We emphasize that this is purely for convenience; our results below are in- 
dependent of this gauge choice. In taking the norm, we use the metric diag(—c 2 , 1, 1, 1) 
i.e. 



|x| 2 = x M x M = — c 2 (x 



°) 2 + £(x<) 2 (5) 



3 We will use the Dirac delta function frequently in this paper. Our presentation maintains the level of 
rigor typically found in the physics literature; in all the calculations below, we always integrate the delta 
function with another well behaved function. However, the reader interested in the issue of products and 



compositions with distributions is referred to 18 1 
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Greek indices run over ... 3, while Latin indices run over the spatial directions 1 ... 3. 



Along the lines of 



191 ]. in this paper we will modify the short distance structure of Maxwell 



electrodynamics by modifying G to 

G(x,y) = S(\x-y\ 2 + d 2 )9(x°-y ). (6) 

In the expression above, the null cone has been replaced by a hyperboloid; this regulates 
the singularities of the theory and leads to a theory analogous to a theory of extended 
particles with characteristic size d. It is curious that for — y 1 ) 2 = 0, equation (J6j) 

forces x° — y° = d. For particles extended in space, one would expect the Green function to 
have support on a spatial interval at zero time. The behavior above is somewhat different 
and we may think of our particles as being extended in time! However, we caution the 
reader against taking this analogy too seriously. The mathematical content of the theory is 
specified by equation ([6]) and this is all we will make reference to below. 

The delta function in ([6]) is written in terms of manifestly Lorentz invariant quantities. 
Second, for a vector satisfying (x — y) 2 = — d 2 , a Lorentz transformation cannot change the 
sign of a; — y°. Hence, the Green function in fl6]) is Lorentz invariant. Consequently, our 
modified theory preserves Lorentz invariance. 

Consider now a point particle of charge q, with worldline specified by a function a M (r) of 
the proper time, r, of the particle. We define the field strength, 

i^ = VM> ( 7 ) 

as usual. 4 Using equations ([3]), ([6]), we find that the field produced by this particle at a 
point x is given by 

F »M) = ~ -—72 J— > ( 8 ) 

47re c (C • a) \ (, • ot J 

where d M = and = are evaluated at the retarded proper time tq, which is defined 
as the solution to 

|x^-^(r )| 2 + rf 2 = 0, (9) 

that has x° > a°(r ). Further, = ^ — a^(r ). For a slow moving particle, the delay 
5 = r — tq ~ -. Equation ([9]) defines a hyperboloid. If this hyperboloid intersects the 



4 Note that this definition ensures that even with the modified Green functions, our theory remains gauge 
invariant 



s 



FIG. 1: Self-Interaction with Modified Propagators 




worldline of the particle at the point 'P' (see Figure [JJ), then the particle when it has 
reached 'P' feels the effect from its past position at 'O'. The equation of motion of a particle 
of mass m is given by 5 

= ^-a u F^ u , (10) 
m 

where F^ u is the net field strength and includes not only the contribution of other particles 
but also a contribution from the past history of the particle itself as we explained above. 
Despite this, the expression for F^ v in (|S]) is free of divergences. Moreover, the self-interaction 
in (|T0|) leads to a finite radiative damping term. 

Second, notice that even for a free particle the equation of motion ([TO]) is in general a 

5 This is the usual Heaviside-Lorentz force law. Our modification of the Green function changes the fields 
produced by a given source but we use the same force law as in Maxwell electrodynamics 
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functional differential equation because it relates ci M (r) on the left to o; m (to) and its deriva- 
tives on the right. To solve this FDE we need the past history of the particle, apart from 
initial data. For FDEs where the delay is bounded above by some So, S < So, the past history 
of the particle specified on an interval of length So (with initial data specified by continuity 
at the endpoint of this interval) suffices to fix a solution 



A. An Example 

Qualitatively, one can understand by a simple linear example how functional differen- 
tial equations help in getting rid of the runaway solutions of the Lorentz-Dirac equation. 
Consider the equation 

m = a w -<)-f® , 

where a, e > are constants. If we expand this equation to first order in e, we find 

(a+l)/'(t) = fr(t)+0(6 2 ). (12) 
If we discard the higher order terms, then (Tl2|) has the solution 

2(l + n)t 

/naive = «1 + K 2 e ™ , (13) 

with K\, k 2 being constants of integration determined by the initial conditions for / and /'. 
Note, that for arbitrary k 2 ^ 0, these solutions grow exponentially in the future. 

On the other hand an exact solution to ffTTj) requires past data over the entire interval [0, ej. 
Substituting, e mt as a test solution in (fTTj) . we find that m must satisfy the characteristic 
equation 

m = ^(e- me - 1). (14) 

It is easy to show that all solutions to ffT4l) have Re(m) < 0. Hence, ffTT]) actually has no 
runaway solutions of the form (TIB"]) . The process of expanding FDEs to obtain ODEs, by 
means of a 'Taylor' series, may lead, as above, to spurious solutions; but this is the process 
that is used to derive the radiation damping term. If e is very small, then locally it is 
indeed true that fit) is well approximated by (TP2l) ; nevertheless the global properties of 
the solution ( 1131) are completely incorrect! This leads us to expect that runaway solutions 
may be avoided by replacing ODEs of the Lorentz-Dirac theory with FDEs, an expectation 
confirmed by the calculation carried out below. 
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B. Many Particles 

To conclude this section, we would like to remark that the formalism discussed here 
works for many particles as well as it does for a single particle. If we have many particles, 
we calculate the field due to each particle by using (jBD. The equation of motion for each 
particle is now given by (fit)]) , where F is the total field strength at a point, obtained by 
summing all the individual field strengths. 

It is notable, that equation (fTUj) is considerably simpler than ([1]) which applies in Maxwell 
electrodynamics. To obtain (JTJ, we need to separate the field strength into two pieces: (a) 
the field strength due to all other particles and (b) the field strength due to the particle 
itself. The second term — when regularized — leads to the second term in (jT^. It is this 
term that is responsible for the runaway solutions that we have discussed. Equation ffTU]) 
involves no such artificial separation. The self action of the particle on itself is regular. 
Since it is this term that leads to runaways in Maxwell electrodynamics and since we do 
not find runaways in our one-particle example below, we will not find runaway solutions in 
a many-body system either. 



III. NUMERICAL SOLUTION FOR A CENTRAL FORCE 

We now turn to a description of how the equations of motion ( JTU1) may be solved nu- 
merically for the important test case of a charged particle moving in an attractive Coulomb 
field. 

We represent the effect of this field by the field strength 

^ r(x) = 4^' ^ Ul W = 0- (15) 

where r = x*Xj. The equation of motion for a particle of charge q moving in this field is now 
given by ffTUj) with 

F» U = F^ + F»f, (16) 
where F self is calculated, as explained in the previous section, with the help of formula (jBD. 
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A. Choice of Code 



To solve the equations (fTOI) . we need a numerical solver for stiff functional differential 
equations. The equations above are stiff for two reasons. First, the time period over which 
we want to solve the equations is very long compared to the delay. Second, the magnitudes of 
the Coulomb force and the radiation reaction differ substantially. Two popular earlier codes 
for FDEs, retard Q and ARCHlf^l] both use the Runge-Kutta coefficients of Dormand and 
Prince. While these coefficients can be used to construct excellent general-purpose solvers 
for ODEs, they are not appropriate for stiff equations. (For an example, of how the dopri 



code fails, see 



HP 



In contrast, the radar code of Guglielmi and Hairer Il5| adapts to FDEs an implicit 
Runge-Kutta (collocation) method based on Radau nodes. The corresponding radau code 



23j is a long-established code which has successfully been applied to various stiff ODEs. 



Another technical issue is that, although we have assumed retarded propagators, the FDE 
( JTOl) is classified as a 'neutral' equation, in the classification of differential equations with 



deviating arguments 



24|, rather than a retarded equation, since the second derivative of 



o: m (to) appears also on the right. Both archi and radar codes may be used also for 
neutral equations, although the radar code is significantly more efficient (by a factor of 
about 10 3 in our experience) in view of its better stiffness properties. The efficient solution 
of such FDEs requires code that can take step sizes much larger than the delay: the radar 
code allows for this possibility. The successful use of the radar code, however, requires 
relatively greater preliminary work to evaluate various Jacobians. 



B. Rewriting the Equations for Numerical Solution 

For the numerical solution it is useful to cast the equations ffTUl) in 3 + 1 form and also 
rewrite them as first order equations. It is easy to see that with the introduction of an 
intermediate three vector u l , and an external time t, we may rewrite the equations as 

m— = qF i0 + q—Fy, 
at 7 
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where 7 is determined by 

7 = a/1 + u.u/c 2 , 
and F is calculated using equations (El), (fTBT) . (fTBT) . 



(18) 



C. Units, Initial Condition and Parameters 

The next task is to choose appropriate units. Since we wish to deal with values that are 
relevant to the hydrogen atom, it is convenient to choose units so that 

c = 300.0, m e = 2.0, -— = 1.5 x 10 3 , q = 0.1, (19) 

47re 

where m e , q are the mass and charge of the electron. This corresponds to a choice of units 
for length, mass, charge and time as 

3.382 x 10" n m, 4.555 x 10~ 31 kg, 1.602 x 10~ 18 C, 
3.384 x 10~ 17 s. 

For the past history, we used a circular orbit with initial velocity appropriate for the first 
Bohr orbit. This is given by v- m = 0.00730c. 

For computational convenience, we will take, unless otherwise specified, d — 1.0 x 10 -3 . 
However, as we explain below, the solution is not very sensitive to the exact value of d. For 
comparison, the classical radius of the electron is do = -^r 2 m c 2 = 0-417 x 10~ 4 in our units. 

The mass that we actually insert into equation ( ITTj) . is taken to be 

q 2 

m = m e — . (20) 

This is because, as explained below, the self-action term, locally contains a part that may 
be thought of as 'renormalizing' the mass back to the value m e . 



IV. RESULTS 



First, let us briefly investigate the solutions to the Lorentz-Dirac equation in a central 
force 



• y TT'COUI I 



mix u = a v F^ + 



6ire Q c 3 



ct „, - a, 



a a. 



(21) 
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This ODE is already stiff because of the large difference in the Coulomb term, the Schott 
term and the radiative damping term. However, it may be solved using the code radau 23| . 
The solution is shown in Figure §2§. It may be seen that the particle shoots off with the a 
velocity that is very close to the speed of light. This is clearly unphysical. 



FIG. 2: Solution to L.D. equation 
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We now turn to the solutions of the equations (TlTj) . Recall that to solve the FDE we 
need to specify the motion of the particle over an appropriate past interval. As stated above, 
we take this past history to correspond to a circular orbit with the parameters of the first 
Bohr orbital. The solution obtained by solving (JTTJ) using the code radar is then shown 
in Figure The damping due to radiation causes the curve above to deviate very slightly 
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FIG. 3: Solution to FDE 
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from the circular orbit which would have been the solution without radiation. This deviation 

is plotted separately in Figure HI Note, that while the particle does spiral in towards the 

center, this is through a complicated elliptic spiral. The average rate of decrease of distance 

towards the center may be measured from the graph in Figure H] and is found to be 

Sr 
~St ' 



-1.13 x 10" 



(22) 



This average value is in good agreement with the naive rate of electron decay, 



Sr 
St 



-Ac 



, , . . — I , which is derived by assuming that the orbit remains circular while losing 
energy at a rate given by the Abraham Lorentz formula in each revolution 

Second, we note that the solution presented above is not very sensitive to the parameter 
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FIG. 4: Deviation from Circular Orbit 
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d introduced in ([6]). In fact, if we change parameters from d 2 = 10~ 6 to d 2 = 10 -7 , while 
choosing m as in (l20l) . we find that the two solutions are almost identical and the deviation 
between them, in the radial coordinate, oscillates with an amplitude of 4 x 10~ n and does 
not increase with time. This may be seen in Figure [5] We also note that this solution is 
in excellent agreement with the backward solution of the Lorentz-Dirac equation discussed 
above. 

Locally, the instantaneous radiation force is described very well by the expression: 

v? 



fi 



rad 



qF^ + q-F^ » jf 



d ( da 1 



d ( d 2 cij 



° l dt I dr ) ° 2 \ dt \ dr 2 



dt c 2 



(23) 
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FIG. 5: Difference between Solutions: d 2 = 10 6 and d 2 = 10 
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where, to an excellent numerical approximation, 

^ 2 2 
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(24) 



87reo dc 2 ' ' 3 

The term c\ effectively renormalizes the mass. This is the reason, we choose m as in equation 
( l20ll . /j rad oscillates as the particle moves in the orbit, and in the Table below we present, for 
different values of d 2 , the root mean square of the deviation (/ rad — / cmp ), which we denote 
by cr, divided by the root mean square of the radiation force (/ rad ) itself, which we denote 



17 



by/rms ' 



d 2 


a 

/rms 


10~ 6 


1.8 x 10- 


-10 


5 x 1(T 7 


1.2 x 10" 


-10 


10- 7 


6.8 x 10- 


-12 


5 x 1(T 8 


7.5 x 10- 


-12 



We see that the fit is always very good and although, initially, the fit becomes slightly better 
as we move to smaller d 2 this gain may be outweighed by numerical errors beyond a certain 
point. 

The stiffness of the equations, and hence the time taken for numerical solution, increases 
as d decreases. Eventually, the numerical solution becomes unstable for a limiting value of 
d, that we find empirically to be d ~ 1.66do, where do is the classical radius of the electron. 

V. CONCLUSIONS 

1. In this paper we discussed the difficulties in obtaining solutions to the Lorentz-Dirac 
equations. We described in Section II CI why these difficulties are especially acute for 
the many-body problem. 

2. Next, we proposed a modification of Maxwell electrodynamics via the modification of 
the Green functions in (161). The modified Green function regularizes the self-action of 
a particle. This leads to the equation of motion (jlOp . A radiation damping term - 
that arises from the self-action - is automatically included in (1101) since involves 
a contribution from the particle itself. 

3. Instaneously, this radiation damping term closely mimics the radiation damping term 
in the Lorentz-Dirac equation. However, globally, the solutions of (fTUj) behave very 
different from solutions to the Lorentz-Dirac equation. We no longer find runaway 
solutions to the equations of motion! A heuristic understanding of how this happens 
is given in Section III Al 

4. We also explained that our modification of Maxwell electrodynamics was (a) Lorentz 
and gauge invariant (b) reduced to the usual Maxwell theory at large distances and 
(c) led to simple equations of motion. 
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5. As an example, we solved our new equations of motion for the important test case of 
a single charged particle moving in the Colulomb field defined by (|15p . Figure [2] shows 
how an attempt to numerically integrate the Lorentz-Dirac equation in this physical 
situation leads to violent runaways. Figure [3] shows the solution of our modified theory 
which is manifestly regular. Figure H] shows that if one looks at this solution closely, 
one finds that the particle is slowly spiralling in, due to radiation losses, as one expects 
on physical grounds. For this simple problem, solutions to the Lorentz-Dirac equation 
may also be obtained using a trick explained in Section II CI The two solutions are in 
excellent agreement. 

6. Having established and tested our technique in this paper, we hope to apply it to other 
concrete physical problems in a forthcoming paper. 



VI. FUTURE DIRECTIONS 



We expect that the technique we have described here will have wide applicability. An 



immediate application o 



physical contexts 



2a 



26|]. 



our technique is to the simulation of relativistic plasmas in astro- 



Another application is to the two body problem that has seen some recent interest. The 
possibility that the classical solutions corresponding to quantum Bohr orbits may have dis- 
tinctive features was discussed in {^J and also conjectured from a very different viewpoint in 
3, J]. These possibilities can now be investigated numerically and may lead to a better un- 
derstanding of the link between differential equations with mixed type deviating arguments 
and quantum mechanics discussed in \^\ . (Recently, a link between pure delay differential 
equations and quantum mechanics was also discussed, from a somewhat different point of 



view, m 



Long ago, the problem of radiative damping was taken up by Wheeler and Feynman in 
the hope that clarity about classical electrodynamics might help to resolve the divergences of 
quantum electrodynamics. In this context, it is noteworthy that the technique of modifying 
the support of the propagators in position space immediately generalizes to a method of 



regularization in quantum field theory, in which context it was originally proposed 19] . The 
regulator here is not merely formal, as in dimensional regularization, and moreover, unlike 
naive cutoffs it preserves Lorentz invariance. This may present some advantages when one is 
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interested in studying the effect of a UV cutoff. In the quantum theory, the unitarity of the 
S-matrix places stringent restrictions on the propagator. This leads to additional difficulties 
that we do not discuss here. 

Finally, in the above discussion we have taken for granted that the interactions involved 
must be retarded, in Q it was shown that the more general equations with mixed-type 
deviating arguments could be used to recover much of the mathematical formalism of quan- 
tum mechanics. In the present context, the use of mixed- type equations may be expected 
to decrease the radiative damping from the full value it has in the retarded case. However, 
solving mixed-type equations involves additional mathematical complexities. We defer this 
discussion to a future paper. 
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Appendices 

APPENDIX A: COMPARISON WITH OTHER METHODS 



In this Ap pen dix, we briefly discuss another popular method, that was proposed recently 
by Rohrlich 29| to solve the problem of runaways. 6 This paper adopts a 'perturbation 
theory' like approach to derive a new equation for a classical point charge. To understand 
this approach, consider a non-relativistic point charge moving under an external force: 



ma 



F ext + 



Q 



da 



(Al) 



6ire c 3 dt 

If we now consider the second term on the right hand side to be a small perturbation on the 
first term, we can rewrite this equation as 

- 2 dF P , 



ma 



F ex t + 



Q 



67re c 3 m dt 



+ ... 



(A2) 



6 After this work was completed, the paper 30| appeared that uses this method to analyze the behavior of 
a charged particle in a Coulomb potential. 
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A covariant version of this equation was derived in [29]. Although the paper initially claimed 



that the resulting equation was exact, this was shown not to be the case in |3l|. Furthermore, 
since the structure of the Lorentz-Dirac equation is tightly dictated by Maxwell electrody- 
namics a modification of the form (IA2I) cannot be performed without modifying the Maxwell 



theory itself. Unfortunately, no such modification was proposed in |29_j. Finally, when rela- 
tivistic effects are included we note that, once again, triple derivatives appear on the right 
hand side of flA2j) which leads once again to third order equations of motion. 



APPENDIX B: NEW FIELD EQUATIONS 

Here we briefly discuss the new field equations that result from the modification of prop- 
agators in equation ([6]). We stress however that this is peripheral to the program that we 
have presented. Given any distribution of particles, we can derive the corresponding field 
strengths using formula (JSj) . We can then calculate the subsequent motion of all particles 
using ( TlOl) . In this formulation of the many-body system in terms of functional differential 
equations, the fields are only an intermediate crutch. 

However, for some formal purposes it may be of interest to analyze the field equations 
themselves. The modified field equations are easily obtained by means of a Fourier transform 
of G. Since these field equations are non-local but yet linear they are more conveniently 
represented in Fourier space and this is what we shall present here. 

To lighten the notation below, we shall put c = 1, in this Appendix. 

Our task is to find the Fourier transform of the modified Green function ([6]). 

GW =/e-»-^ + d *) W (Bl) 

Formally, this integral is identical to the Fourier transform — from momentum space to 
position space — of a massive propagator. Using the standard techniques, developed for 



this (see, for example [32|), we find 



(2tt) 3 4vr ' 1 8nyf-p 



lG{, ' ] 1 s(p°)6( P 2 ) + idv\-%£± \Y X (dy^) - isip )^ U^Ttf 

I R77-1 / — n 2 L V / V 

K\ {m\fp 



6{p 

47T 2 v 



;b2) 
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where J and Y are Bessel Functions of the first and second kind, K is a modified Bessel 
function {33I and s(p°) is the sign of p°. All terms above should be understood to be 
distributions and V denotes the principal value. Note that the expression (1B2[) is Lorentz 
invariant. This is because, for p 2 < 0, we cannot change s(p°) by means of a Lorentz 
transformation; all other terms in (1B2[) are written in terms of explicitly Lorentz invariant 
quantities. 

Let us consider this expression near p 2 = 0. In that limit, flB2j) may be written as: 



-iG(p) 1 2 
— = -s{p )8(p ) - V 



Att 2 p 2 8 



id 2 j d^/\p 2 \ 



71- 



V 



-d 2 

16n 



s(p°)o(-p 2 ) + o (-/FfiB b 2 



(B3) 



In the physics literature, the first two terms above are often combined 
This also tells us that, in the limit d 0, 



4ir 2 (p 2 — iep° 



(B4) 



lim G(p) 



2tt 



(B5) 



which is familiar as the usual expression for the retarded Green function in momentum 
space. The new field equations are now simply 



(B6) 



where A(p) and j(p) are the Fourier transforms of the vector potential and the current 
density. 

A» = j A^e-^d^x, j» = j j^e-^d^x (B7) 

We note, that in the limit that d — > 0, as a consequence of (1B5I) . these equations reduce to 
the equations of Maxwell electrodynamics 



p 2 A(p) = j(p) 



(B8) 
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